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ABSTRACT 


One method of reducing truncation error in numerical 
prediction is the use of more accurate finite-difference 
approximations. A relatively simple barotropic model, 
Similar to the one currently in use by the Fleet Numerical 
Weather Facility (FNWF), Monterey, California, is employed 
as a basic program in testing several higher order finite- 
differencing schemes. Forecasts were computed to 24 and 
48 hours with modifications of the basic model, and 
comparisons were made with the analyses at verifying time. 
It is the intent of this study to determine the effects of 
some higher order finite-difference approximations in this 
numerical prediction model. 

The author expresses his appreciation to Professor 
George J. Haltiner for his suggestions, guidance, and 
patience. Appreciation is also expressed to the personnel 


of FNWF for providing data and programming assistance. 
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i imeroduction:. 

In any numerical method of weather prediction, 
regardless of the simplicity or complexity of the forecast 
model, one source of error which cannot be entirely 
eliminated is the so-called truncation error. This error 
is mathematical in nature and stems from the use of finite- 
difference approximations to the derivatives of functions 
such as those encountered in meteorology. Thompson [7] 
discusses the problem of truncation error and outlines 
two methods of reducing the magnitude of the error. The 
first method involves decreasing the grid interval thereby 
increasing the number of grid points for a given area. 

This procedure is not considered in this study. The second 
method utilizes the "quasi-analytic" property of the 
functions dealt with in numerical prediction which implies 
that the functions and their derivatives vary smoothly. 

If in reality these functions are analytic, uniformly 
convergent Taylor series can be written for any such 
functions. The inclusion of higher order terms of the 
Taylor expansion results in more accurate finite-difference 
approximations. 

It is the purpose of this study to determine the 
effects of some higher order finite-difference approximations 
in numerical weather prediction with the intention of 


obtaining a more accurate prognostication. 





is The Basic Model 

The basic model used in this study is a simple model 
described by Haltiner [4] and similar to the 500-mb 
barotropic-type forecast currently in operational use by 
the Fleet Numerical Weather Facility (FNWF), Monterey, 
California [2]. This model was selected because of its 
relative simplicity, availability of data in a convenient 
form, and suitability in programming. It was supposed that 
a simple model would lend itself best to this type of study. 
In this model, three empirical terms are incorporated; 
(a) a "Helmholtz" term, (b) a "latitudinal" term, and 
(c) a "temperature" term. Although these terms are 
reported to be relatively small, they are included because 
they have been found to improve the results. 

The equation programmed for the basic model is of the 


form: 


(yp) 22 =~CJ(zZ,2+H). - 


Equation (1) is solved numerically by a sequential 
relaxation method. The input parameters are the 500-mb 
analysis and the thickness of the layer between 1000 mb and 
500 mb, with only the values north of 2N entering into the 
computations. A forward time step is used initially 
followed by centered-differences at one-hour intervals. At 


six-hour intervals, slight smoothing is applied and the 





ellipticity criterion imposed, after which a forward time 
step is again initiated and the process repeated. 

The finite-difference approximations programmed in 
the basic forecast scheme can be found in most texts on 
numerical methods. Collatz 11] provides a comprehensive 
list of these formulas. Formulas used in the basic program 


are the 5-point approximation to the Laplacian 
zoe A (Zz, +2,+25+Zy-4Z, ) 
ee l 2 5 4 o 


and the 5-point approximation to the Jacobian 


J (4) = in| (t- bY Y,-%)-(%- H- i) 


where the numerical subscripts refer to the grid position 


with the 5-point grid defined as 
“ @3 
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The time stepping formulas consist of the forward-difference 


eeawar ee (22) at 


(4) 
and the centered-difference 
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where the subscripts refer to the time at which the 
quantities are computed. 

To reduce the truncation error in the basic forecast 
program many higher order finite-difference formulas are 
available. Only a few of these formulas are considered 
in this investigation. To identify the modifications to 
the basic program, letter designations were given to the 


changes which are described in the following section. 


or Modification B. 

The first modification to the basic program consists 
of substituting for equation (2) a 9-point approximation 
to the Laplacian 

UZ, = gya[ 602. +16(2,+2,+2,+2)- (Bt ZB]. 503 
Equation (6a) can also be written symbolically as suggested 
by Forsythe and Wasow [Bales 

=i 
/6 
Y ga -| /6 ~ fone 
zi (6b) 
The remainder, representing the residual between the 


finite-difference approximation and the actual value, 


consists of sixth-order derivatives and higher. 


ive Modification C. 


The same grid as presented in section 3 can be used to 
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derive a 9-point approximation to the Jacobian applying 
the formula for a 5-point centered-difference along each 
dimension of the grid. Modification C includes the 


replacement of equation (3) by 


Tb) = | (3-64 +86,-$,)(¥, -8h,+82-$,) 
a ($,- 8h +8%,- & )(, - ay ap ~ .) | 


2) 


where the subscripts locate the positions on the grid as 


follows: : 
e 
e~ 
Ss ' 0 2. 6 
€ @ © % ° 
ett 
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The remainder contains the product of third-order derivatives 
as the leading term. In both modifications B and C no other 
changes are made to the basic program, and no additional 
computer storage space of any significance is required. 

By the inclusion of both 9-point approximations, equations 
(6) and (7), a composite program is obtained. This dual 


program is designated as modification X. 


Ds Modification D. 
In the basic program, equations (4) and (5) are used 
for the first and second hourly time steps respectively. 


For this modification, subsequent time steps are performed 
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by utilizing an approximation comprised of terms through the 


third derivative as shown: 


3 
Zine eee 2 (22) ace & (28) (ety. 
(8) 
By the method of divided differences described by Jennings 
eB a formula for equation (8) in terms of 2 and Sécan be 


derived: 


ore RS ($= ),4 = Be. tO 2 ae - z age 
(9) 

Equation (9) has the same form as the formula listed in 
Milne [6] for a net of four equally-spaced points. The 
remainder can be shown to contain terms higher than the 
third derivative. 

Note that equation (9) contains one additional past 
point. In programming this modification, one additional 


field is necessary in the computer memory. 


oe Modification ©. 

Iterations after the first two time steps are 
accomplished in this case by an alternative method of 
including the third-order term in equation (&). By 


considering 
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and using a one-sided difference for the second derivative 


in terms of the first derivative, equation (8) becomes 


i ‘ (SF), - ity: eo GP)..4 | 

(10) 
Equation (10) contains three height-change tendencies and 
consequently greater computer storage space is required. 
Due to limited computer storage space and other considerations 
to be pointed out in section 7, both modifications D and E 


necessitated considerable revision of the program. 


y Procedures. 

All computations were performed with the CDC 1604 
all purpose computer using a 63 x 63 square grid with a mesh 
length of 381 km. The grid is superimposed on a polar 
stereographic projection of the northern hemisphere centered 
at the North Pole and true at 60N. 

Imput data in the form of 500-mb and 1000-mb analyses 
on magnetic tape were provided by FNWF. The period from 
O00Z 15 March 1964 through 00Z 20 March 1964 was used in 
this study. This period was selected because the synoptic 
situation showed a block positioned over northern Europe, 
a series of waves moving rapidly across the Pacific area, 
and systems developing and weakening across the North 
American continent. This depicted a variety of conditions 


for which the programs could be tested. 





The prognosis was carried out to 24 and 48 hours for 
all programs with the exception of modification E. 
Difficulty was encountered with modification D because of 
overflow after 36 hours. By a small amount of smoothing 
at each time step the program was able to operate to no 
less than 48 hours. Modification E overflowed at 30 hours. 
smoothing the height field or the tendency field was in- 
effective in controlling the instability, hence modification 
Kk was terminated at 24 hours. 

When each series of computations was completed, the 
products were compared at each grid point with the analysis 
at the appropriate verifying times. The pillow and root 
mean square error were computed for each prognosis. The 


pillow and root mean square error are denoted respectively: 


and 


Rio e 


if 





where n is the number of grid points considered. For this 


comparison all grid points were used, hence n=3969. 


S. Results and Discussion. 
The results of the comparisons are given in Table I. 
It should be noted that the results are slightly biased in 


favor of the basic program. This is due to the use of the 
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TABLE I 


Pillow and RMSE of verifying analysis and prognosis 


2h—-hr 4é~-hr 
Prognosis Initial Time oie R/RMSE 
Basic 15/002 6. yom 7 7.9/50.3 
16/002 25 yee 3.7/0 
17/00Z 2.1/35.4 0.3/54.9 
18/002 105) 20)45 LatS/ipoe > 
Mod B 15/00Z (20/4 Ia 34) SOT 
16/002 37 GOO Be Gee 
Mod C 15/002 ey Sly) 9.4/49.1 
16/00Z 3.3/30.5 Bel ee 7 
17/002 oa ee Fre ao 
18/002 0.3/29.9 ae poe 
Mod D 15/002 16.8/48.8 44.2/97.6 
16/00Z T2527 pore 39.6/86.0 
Mod E 15/002 G56) 02m = 
16/002 4s 3 (bere - 
Mod X 15/002 73/7 50.0 8.8/68.0 
16/002 ory 7 0 5.5/74.0 


basic prognosis in the first approximation for producing 


the verifying analysis. Maps containing the 500-mb height 
field were printed as a grid for each set of computations 


and contour charts plotted with the CDC 160 and graph plotter. 
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The prognoses from 00Z 16 March 1964 for each program, 
typifying the products obtained, are given in Appendix I. 
Analyses of the O00Z charts covering the period concerned are 
also presented. 

By examination of the products of the various programs 
and comparison with the analyses at verifying time, the 
following observations were made: 

(a) All programs show reasonable contour patterns 
with the exception of modification D. It is obvious from 
charts 7 and 8 that spurious perturbations with a wave- 
length on the order of three grid-intervals developed and 
amplified with time. These waves appear to be caused by 
computational instability originating from the use of the 
finite-difference formula. 

(bo) The major polar low was too deep in all programs, 
ranging from 20-60 m at 24 hours to 105-160 m at 48 hours. 
On the other hand, the heights of quasi-stationary systems 
tended to be too large by an average of 12 m for highs and 
70 m for lows at 24 hours, and 43 m for highs and 71 m for 
lows at 48 hours. The magnitude of the differences tended 
to contribute to the relative accuracy of the individual 
programs. 

(c) The speed of moving systems was in general slower 
than the verifying analysis indicated. Interestingly, move- 
ments in modifications B and X could hardly be detected in 


most cases. Recall that modifications B and X both contain 
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the 9-point approximation to the Laplacian. This change 

was suspect to possible programming error or incompatibility 
of the modification and the basic program. A check of the 
subroutine for the 9-point Laplacian gave reasonable 
calculations. Consequently, the latter reason is indicated 
as the cause of the results, especially in view of the 
tendency for spurious anti-cyclogenesis to creep in. Anti- 
cyclogenesis was controlled in the other programs. 

(d) As expected for barotropic models, all programs 
suffered when predicting the intensification of developing 
systems. 

Of all the modifications tested, the only one which 
compares favorably with the basic program is modification 
C which utilizes the 9-point Jacobian. There is strong 
evidence that the use of this subroutine can improve the 
basic program. Due to the limited computer time available, 
the length of each run (see table II), and other difficulties, 
a relatively small sample was obtained. Because the 9-point 
Jacobian program showed more promise, additional data were 
collected for modification C. Similar data for the basic 
program were also required to determine if the two methods 
were comparable. The RMSE for modification C averaged 
slightly less than the basic program, though the pillow 
was somewhat higher. There was no significant difference in 


running time. 


er 





TABLE II 


Computer running time for the prognoses 


Prognosis 24-hr Pro 48-hr roe 
fe Oe minutes 


Basic i 24 
Mod B oe ae, 
Mod C LZ 2h 
Mod D 2k 47 
Mod E ay, - 
Mod X 16 32 
9. Conclusions. 


From this investigation, it can be said that although 
truncation error can be reduced theoretically through the 
use of higher order finite-differences, the incorporation 
of these approximations into a numerical prognostication 
scheme can lead to undesirable effects in some cases. 
Computational instability and other difficulties encountered 
in this study point out the implicit complexity of making 
a seemingly simple change to a forecast method. No 
modifications other than the substitution of higher order 
finite-difference formulas were used, except for the small 
amount of additional smoothing previously mentioned. By 
adding empirical corrections to the model, it is conceivable 


that the forecasts could be further improved, but this was 


eZ 





not the primary objective of this paper. 

The limited results of this study provide additional 
avenues of research. One-sided advection schemes and 
numerous other stencils for the finite-difference approxi- 


mations are suggested for further examination. 
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APPENDIX I 


The following charts were produced during this investigation. 
They are typical of the prognoses obtained by the programs 
studied. Also included are the analyses of the charts for 
the period covered. 
Charts 1-2 24-hr and 48-hr prognoses, OO0Z 
16 March 1964 Basic program 


Charts 3-4 24-hr and 48-hr prognoses, OOZ 
16 March 1964 Modification B 


Charts 5-6 24-hr and 48-hr prognoses, O00Z 
16 March 1964 Modification C 


Charts 7-8 24-hr and 48-hr prognoses, OOZ 
16 March 1964 Modification D 


Chart 9 24-hr prognoses, OOZ 16 March 
1964 Modification E 


Charts 10-11 24-hr and 48-hr prognoses, O00Z 
16 March 1964 Modification X 


Charts 12-17 Analyses, O00Z 15 March 1964 
through OOZ 20 March 1964 
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